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Splitting integrators for the BCS equations of superconductivity 
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The BCS equations are the centerpiece of the microscopic description of superconductivity. Their 
space discretization yields a system of coupled ordinary differential equations. In this work, we 
come up with fast time evolution schemes based on a splitting approach. One of the schemes only 
requires basic operations. For the physically important case of the BCS equations for a contact 
interaction potential, the computational cost of the schemes grows only linearly with the dimension 
of the space discretization. Their accuracy is demonstrated in extensive numerical experiments. 
These experiments also show that the physical energy of the system is preserved up to very small 
errors. 
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1. INTRODUCTION 

In this work, we consider the time-dependent BCS 
equations, often also referred to as Bogolubov-de-Gennes 
(BdG) equations. These equations, named after Bardeen, 
Cooper and Schrieffer , see [T|, are are the basis of the 
microscopic description of superconductivity in metals. 
They are coupled partial differential equations which de¬ 
scribe the evolution of the particle density 7 and the 
Cooper pair density a of a fermionic system (see also [2] 
for a detailed introduction). Although the BCS/BdG 
equations are a fundamental part of condensed matter 
physics, their numerical treatment has not been paid at¬ 
tention to so far. Thus, in this work, we present two re¬ 
liable and efficient integration algorithms for these equa¬ 
tions based on a splitting approach. 

The evolution equations for the Cooper pair density 
resembles the linear Schrodinger equation for quantum 
dynamical systems. One important aspect of these sys¬ 
tems is that, after a space discretization, the right hand 
side of the resulting ordinary differential equations has 
a very large Lipschitz constant caused by the Laplacian 
in the kinetic part. As a consequence, standard explicit 
integration schemes, such as the ones presented in [3], 
although very popular in computational physics, are of 
no use in quantum mechanical applications. Therefore, 
the treatment of such quantum dynamical systems has 
been of huge interest in the numerical analysis commu¬ 
nity for many decades, see, e.g., [H Chapter II. 1]. Vari¬ 
ous evolution schemes for the linear Schrodinger equation 
in varying settings have been proposed over the years, 
see, e.g., [MID]- Nonlinear Schrodinger equations such 
as the Gross-Pitaevskii equation and equations arising 
from the Hartree and Hartree-Fock approximation of the 
quantum state have also been devoted attention to, see, 
e.g., [Mill! and [HI US], respectively. 

All theses methods have in common that the partial 
differential equations are first discretized in space. This 
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means that the system is restricted to a suitable subspace 
spanned by a finite number of basis functions. Here, we 
do this with the help of a Fourier collocation method 
which is the straightforward approach for the problem 
we look at, see, e.g., [U Chapter III. 1]. This yields a 
system of coupled ordinary differential equations, on the 
solution of which we focus in the present work. 

What, in many applications, turned out to be the 
most promising tool for the solution of the space dis¬ 
cretized system was the splitting of the equations under 
consideration into some subproblems, each of which can 
be solved more easily than the system of equations as a 
whole. This idea was first employed for advection equa¬ 
tions in US and nu. In the realm of quantum dynamics, 
it was applied for the first time in m where the lin¬ 
ear Hamiltonian was split into a kinetic and a potential 
part. The respective solutions were then concatenated 
in a suitable way in order to obtain a reliable integra¬ 
tion method. Here, we use this ansatz to introduce two 
schemes for the evolution of the space discretized BCS 
equations. The coupling terms depend on the convolu¬ 
tion of the particle density with the Cooper pair den¬ 
sity. We use the fast Fourier transform (FFT) to swiftly 
compute these terms. As a consequence, the CPU effort 
per time step of our schemes grows only mildly with the 
number of basis functions. This is very important since 
in most physical applications the BCS system requires a 
discretization space of very high dimension. 

For the first scheme, we exploit that the eigenvalues of 
the density operator, which are functions of the particle 
density and the Cooper pair density, are conserved along 
exact solutions to the BCS equations. Hence, we can ex¬ 
press the particle density as a function of the Cooper pair 
density. We end up with a decoupled nonlinear system 
for the evolution of the Cooper pair density a. The thus 
obtained equations are split into a linear part, which can 
be solved exactly, and into a nonlinear part, the flow of 
which can be approximated by some standard numeri¬ 
cal scheme. In the rest of this work, we will refer to 
the resulting integrator as BCSInt. It is very accurate 
and preserves the physically interesting eigenvalues of the 
density operator by construction. The integrator has al- 
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ready been employed in a numerical study of the physical 
behavior of the BCS equations in m- 

For the second integrator, we do not decouple the sys¬ 
tem at all. Instead, thanks to the system’s particular 
structure, we can aptly split it into three subproblems for 
which the flows can be calculated very efficiently. These 
calculations require only basic operations. Recombining 
the thus obtained flows in a suitable way results in a 
very accurate and efficient scheme which conserves the 
system’s constants of motion, such as the energy, up to 
very small errors. In the following, we will denote the new 
scheme by SplitBCS. In the physically important case of 
a contact interaction, i.e., when the potential is given by 
a delta function, the flows of the subproblems can all be 
calculated exactly with an effort linear in the number of 
basis functions. 

We demonstrate our integrators’ favorable behavior 
with the help of numerical experiments and numerical 
comparisons to standard integration schemes. We men¬ 
tion that an error estimate similar to the one in [Z 01 for 
splitting schemes applied to the Gross-Pitaevskii equa¬ 
tions is expected to hold for SplitBCS. However, such an 
analysis is out of the scope of this work. 

Our presentation is organized as follows: We start with 
a short introduction to the BCS equations in Section [2] 
Afterwards, we explain the Fourier collocation for the 
partial differential equations in Section [3] We first in¬ 
troduce our splitting scheme BCSInt for the decoupled 
nonlinear system in Section [5j Then, we present our fast 
integration scheme SplitBCS for the coupled system in 
Section [ 6 ] This is followed by numerical tests in Sec¬ 
tion [TJ Finally, we summarize our results in Section [ 8 ] 


momentum space representations 

7 (t,P) = ^~ f x)e lpx dx, (3) 

a(t,p) = [ a(t,x)e lpx dx. (4) 

J r 

In this basis, the equations can be written in the compact, 
self-consistent form 


ir(t,p) = [H r ^ tP) ,T(t,p)] , peR, 
see, e.g., { 22 ]- T(i,p) is the 2 x 2 -matrix 

r(t ’ p> = (sm i-wrij 

and Hp(t, p ) is the Hamiltonian 

()= f P 2 Z P 2[V*a](t,p) 
(t, p ) \^2[U*a](f,p) P~P 2 

Here, * denotes the convolution of V with a(t,p). 


(5) 

( 6 ) 

(7) 


2.1. Superconductivity 

It can be shown, see e.g. m, that the free energy 
functional 

F T (T(t))= j (p 2 - p)^(t,p)dp + [ \a(t,x)\ 2 V(x)dx 

JR JR 

+ [ Tr C 2 (r(p)logr(p))dp (8) 

J R 


2. THE BCS EQUATIONS 


A superconducting translation invariant system in one 
spatial dimension is characterized by the particle density 
7 :MxRh)R which describes the probability at time 
t of finding a particle at position x and the Cooper pair 
density a:IxIqC which gives the probability at 
time t of having a Cooper pair of electrons at distance 
x. For a given particle interaction V, the evolution of a 
and 7 is governed by the BCS equations, sometimes also 
called Bogolubov-De-Gennes equations, 


x) = —2 / V(y)lm a(t,x - y)a(t,y) dy, (1) 


i a(t, x) = 2 ( — -j -2 — p + V(x) ) a(t, x) (2) 


-4 / ”/(t,x - y)V(y)a(t,y)dy, 
J R 


with y denoting the chemical potential of the physical 
system and ‘ = d /at. Conventionally, the BCS equations 
are given in terms of the Fourier transforms, i.e., the 


is conserved along solutions of the evolution equations ^ 
for any given temperature of the system T. If, for a 
given temperature T, the minimizer T of J~r has a non 
vanishing Cooper pair density a, then the system is said 
to be in a superconducting state. 


2.2. The discrete BCS equations 

In order to render the system computationally palpa¬ 
ble, one restricts it to a domain D = [0,T27t], L £ N, 
and assumes periodic boundary conditions. In most ap¬ 
plications, L is a large integer as the extension of the 
system is considered to be huge compared to the micro¬ 
scopic scale which here is 0(1). On the finite domain D , 
the momenta consist of the discrete set k £ x liJL. The 
momentum space representations of a and 7 are given by 

-1 pL2,7T 

%(t) = J o 7 (t, x)e ikx dx, (9) 

-1 pL2n 

= a(t, x)e lkx dx. ( 10 ) 

-^ 7r Jo 
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In terms of these representations, the BCS equations read 


3.1. System of ordinary differential equations 


iT*(t)= [H rk(t) ,T k (t)] , ( 11 ) 

where the convolution appearing in the Hamiltonian is 
now to be understood as 

= ( 12 ) 

iez 

The first step to a numerical solution is to introduce a 
finite basis. This process is called space discretization. 


3. SPACE DISCRETIZATION 

As the BCS equations are given in their momentum 
space representation anyway, it is most convenient to use 
the so-called Fourier collocation. This means that for a 
fixed number K G N, a L27r-periodic function f(x) = 
E je z fti) elk/LX is approximated by 

K-l 

f K (x)= E (13) 

k=-f 

where the coefficients fj^ are obtained by the discrete 
Fourier transform of the values fj = f ( L2t /k ■ j), j = 
— k /2 ,..., K j‘i — 1. From numerical analysis, cf. j¥| Chap¬ 
ter III. 1], it is known that for an s-times differentiable 
function /, the bound 

\\f(x)-f K (x)\\<CK- s \\^l\\ (14) 

holds for some constant C independent of the number of 
basis functions I\. 

Mathematically speaking, we work on the subspace 
spanned by the first K eigenfunctions of the Laplacian 
on [0,L27 t]. The approximation of the particle density 
on this subspace is given by 


o ± 


7 K (t,x)= E 7(15) 

k=-f 

and the approximation of the Cooper pair density reads 

i£_i 

E (t)e l i x . (16) 


, K / 




Inserting this approximations into the infinite dimen¬ 
sional BCS equations (111 yields a finite dimensional sys¬ 
tem of ordinary differential equations (ODEs). 


The system of ordinary differential equations we end 
up with after applying the Fourier collocation is given by 

m-W = a k (t)(v * a) - a k (t) (v * cx) , (17) 

idfe(f) = 2 ~ t^j a k(t ) - (2y k {t) - 1) (v * , 

(18) 

K , I< i 

- <k< -1, 

2 “ “ 2 

where, for the sake of readability, we have replaced *y K 
and a K by 7 and a, respectively. 


3.2. System for a contact interaction 


For a contact interaction V(x ) = ~aS(x), a > 0, which 
is the most popular interaction model in physics, we have 


V(k) = 


a 

2 LE 


K 

T 


< k < 



(19) 


Hence, the convolution term in the self-consistent Hamil¬ 
tonian on the K dimensional subspace is given by 


(V*a K ) (t) = --E 
V / k K 2Ln 


2 x 

E 


af{t). 


( 20 ) 


With this relation, the equations of motion become 


!f/2-1 K/ 2-1 \ 

ijk(t) = ( a k(t) E - ak{t) e > 

j--K / 2 j=-K / 2 J 

( 21 ) 

«/a-l 


idfe(t) = 2 (^2 — a k (t) + -j— E ( 2 7fc(i) - !) > 

V ' j =- K /2 

( 22 ) 

K , K 

—— < k < —— 1. 

2 ~ ~ 2 

With 

p k (t) := Rea k (t), (23) 

q k (t) :=lm.a k (t), (24) 

we can rewrite the equation of motion for 7 k (t) as 


7 k(t)= I q k (t ) E PjiQ-Pkit) E I ■ 


Ff/2-1 


’C/l-l 


Ltt 


j = - K /2 


j=~ K /2 


(25) 


From this expression we can see very easily that 7 k (t) is 
a real quantity whenever 7 t( 0 ) is so. As 7 represents the 
physical particle density, which is real by definition, we 
can safely assume 7 k (t) to be real in the following. 
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3.3. Constants of motion 


4. CALCULATING THE CONVOLUTION 
TERMS 


For later use we mention that the coupled sys¬ 
tem (17 ),(18 1 possesses some important constants of mo¬ 
tion: 


It can readily be seen that the matrix -ffr(t) hr the 
BCS equations (111 is self-adjoint. Together with 
the commutator structure of the equations of mo¬ 
tion ( 111 , this implies that the evolution of T(f) 
is unitary. Consequently, its eigenvalues are pre¬ 
served along the evolution. A little bit of algebra 
shows that these eigenvalues are given by 




KWI 2 . 


(26) 


We denote by J- the Fourier transform of a vector of 
length K = 2 N , N £ N, and by J 7 ” 1 its inverse. With 
the help of the fast Fourier transform (FFT) algorithms, 
these operations can be calculated efficiently in 0(N ■ I\) 
operations, see, e.g., m Chapter 12]. 

Furthermore, the convolution of two K dimensional 
vectors a and b can be computed by 

a * b = T - 1 ((Tcl) ■ (Tb)) , (29) 

with • denoting pointwise multiplication. Taking this 
into account, we can efficiently calculate the convolution 
terms as outlined in Fig. [T] There, we have defined 

V r .= V{^n/ K .j), j = — k /2, ..., K /2 — 1. (30) 


• The discretized analog of the free energy func¬ 
tional ([ 8 ]) in the case of an interaction potential 
is given by 


F K h(t),a(t)) 



- M 7 k(t) 


1 u 71 ' 

+ 2Wo ^)H^)I 2 ’ 

K/2-1 

+ T J2 [A+ log(A+) -F log(A-)], 

k=— K f 2 

(27) 


and can be shown to be preserved, too. 


The algorithm only takes 0(N ■ K) operations. When 


Algorithm 1: calc convolution 


conv = FFT(a) 
for j = — K fi to K fi — 1 do 
j^convj = conv, • Vj. 

conv = invFFT(conv) 


Figure 1: Sketch of the algorithm calc_convolution, which 

uses the FFT and its inverse to efficiently calculate the con¬ 
volution between a and V. 


considering a system with contact interaction, i.e., when 
integrating the evolution equations (251,(221, the convo¬ 
lution terms are just a sum over the entries of the vector 
a. Hence, the CPU effort in this case is only O(K). 


3.4. Numerical notation 


5. NONLINEAR SPLITTING INTEGRATOR 


From a numerical point of view, the coupled sys¬ 
tem (]l7|),(fT8j), when supplemented by some initial data, 
represents an initial value problem 


BCSInt, the integrator we present in this Section, is 
based on the conservation of the eigenvalues of F(t). 
These eigenvalues being conserved, the following equality 
holds 


^ =/(y(t)), 
y(0) = y 0 , 


(28) 




+ M )\ 2 



0 2 + Mo)| 2 . 

(31) 


with y £ C 2K . Formally, the aim of this paper is to find 
a numerical approximation to the exact flow of such an 
initial value problem. For this, we denote a time step by 
r and the flow over such a time, i.e., the smooth map 
between y (t) and y (t + r), by <h Ti /(y(f)). Its numerical 
approximation will be denoted by 4)™ 1 . 

Both of the numerical flows we present in this work rely 
on the fast calculation of the convolutions appearing at 
the right hand side of the equations of motion (|l7|),(|T8j). 
Let us turn towards this now. 


With the help of this relation, we can eliminate Jk(t) in 
the equations of motion for ak(t) as we show now. 


5.1. Decoupled system 


Solving Eq. (31) for 7 k yields 


7 k(t) = VHk) - l«fc(f)| 2 , 


(32) 
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with the auxiliary function 


h(k):= (7*(0)-§) +K(0)| 2 . (33) 


The sign in relation (32 ) can usually be inferred from 


physical information. In our study m, for example, the 
initial values had to be such that 7 fc( 0 ) was greater than 
!/2 for ^ > fc 2 /i 2 and less than or equal to 1/2 for /r < 
k 2 / L \ 


Inserting the just-derived expression (321 for 7 *. (t) into 


the equations of motion for a k (t ), we get the nonlinear 
system 


ior fc (*) = 2 ( - n ) a k (t) 


T/ 2-1 


±—^h(k)^\a k (t)\z Y, a ‘O')- (3 4 ) 

j = - K /2 


K K 

- <k< - 1 . 

2 “ “ 2 


Having decoupled the system, we can now turn towards 
its time evolution. 


5.2. BCSInt 


The nonlinear system (34), together with some suitable 
initial data, gives an initial value problem 


^ = Km, 

a( 0 ) = d? 0 , 


(35) 


for 

a=(a_K/ 2 (t ) ... aj f / 2 _ 1 (t)) T e C K . (36) 

The right hand side of the differential equation can be 
written as the sum of two terms, 


f(a) = Aa + (37) 


where f± represents the nonlinear term and where A is 
the matrix 


A = diag 



-,2 



(38) 


This linear part resembles the kinetic part of the linear 
Schrodinger equation. Its flow 4> TjJ 4 can be calculated 
exactly as 


$r,A(a) = diag 



4 L 2 



T 


—i2 

e 


( 


{K — 2) 2 
4 L 2 



<5. 


( 39 ) 


With regard to fi, it has a much smaller Lipschitz con¬ 
stant than the complete right hand side /, wherefore 
can be approximated by some standard integration 
scheme. We than follow the idea of ESI and set 

*"7(S(0)) = (*r/ a ,A o KT° («(«))• (40) 

Applying this operation successively yields an approxi¬ 
mation to the exact solution at times t = nr, n = 1 , 2 ,... . 
Its error decreases quadratically as a function of the step 
size r as long as T" 1 ™ is a second-or-higher order approx¬ 
imation to see, e.g. [Z3 Chapter II.5]. 

Just as every exact flow, satisfies 

®t,A ° ®s,A = ®t+s,A- (41) 

Hence, when applying many time steps of the numerical 
scheme in a row, one can combine the last sub-step of the 
previous step with the first sub-step of the next step, thus 
saving computational costs. We illustrate the resulting 
procedure in Fig. [ 2 ] 


Algorithm 2: BCSInt 


a = 4>7 2 ,A(ao) 

for n = 0 to N do 

a = («)■ 
a = 4> T ,A(a). 


Figure 2: Sketch of our algorithm BCSInt which for a given 
initial value ao and a given step size r approximates a(Nr) = 
4 > jvt ,/(« o )- 


Concerning 'h"’™, in the study [15] it has been calcu¬ 
lated via the fifth order explicit Cash-Karp Runge-Kutta 
scheme proposed in [3], In the experiment Section [7] be¬ 
low, we will also test the second order explicit midpoint 
rule. In this case, d*"™ is calculated as outlined in Fig. [3l 


5.3. Number of operations 

In order to analyze BCSInt’s efficiency, we count the 
number of real operations which are executed per call of 
our implementations, which, to the best of our knowl¬ 
edge, have been implemented in the most efficient way 
possible. We do not weight the costs of different op¬ 
erations, i.e., the square root in calc/i, cf. Fig. [3j also 
counts as a single operation. The number of operations 
as a function of the number of basis functions K for the 
various sub-algorithms and BCSInt as a whole are listed 
in Tab. |T] We mention that, if we substitute the fifth or¬ 
der Cash-Karp scheme for the explicit midpoint rule in 
calcd*^, the number of operations for calcd)^ increases 
to 6 x calc_convolution + 38A'. 

Let us now introduce our second integration scheme. 
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Algorithm 3: calcTy, Algorithm 4: calc/i 


Y =calc/i(a) 

for k = — K /i to K /i — 1 do 
Y(k) = 

a k (t) + r/2Y(fc) 

Y =calc/i(Y) 

for k = — K 1 * 2 . to K /i — 1 do 
a k (t) = a k {t ) + t Y (k) 


c = calc_convolution(V, a) 

for k = — k /2 to K /2 — 1 do 
d = \Jh(k) - |afc(t )| 2 
fi (fc) = —ic • d 


We will now show that we can efficiently calculate the 
flows for all three subproblems. The calculation of <f> T 4 is 
nothing other than < 1 > T .4 acting on a with 7 held constant. 
We thus, in fact, only have to consider the other two 
subproblems. 


6.1. Calculating <f> Ti9 


For the subsystem 


Figure 3: The left panel shows the algorithm which for a 
given value a(nr) and a given time step r calculates <5((n + 

1)t) = <3>" u ?‘(d(m")) with the explicit midpoint rule. The , . , 

right panel shows the algorithm which for a given value d the rl S ht hand Slde does not de P end on the quantity to 
calculates /i(d) be ev °lved. Therefore, the solution of the initial value 

problem (|48|) at time t is trivially given by 


=9(a(t)), 
.7(0) =7o, 


(48) 


Algorithm 

^Operations per call 

Calculation of 3 > t ,a 
C alc/i 
calc*i> 

BCSInt 

BCSInt for contact interaction 

14 • K + 14 

1 x calc_convolution + 12 • K + 20 

2 x calc/i + 8- A + 9 = 2x calc_convolution + 32 • K + 49 

2 x calc_convolution + 46 • K + 63 

58 • K + 63 


Table I: The required number of operations per step as a 
function of the dimension of the ODE system (34 1 for the 
sub-algorithms of BCSInt and for BCSInt itself. 


l{t) = $t, 9 ( 7 ( 0 )) = 7 ( 0 ) + t ■ g{a{ 0)). (49) 


Bearing in mind the reformulation (251, we calculate a 
step of 4> r g with the algorithm illustrated in Fig. [4] 
Please note that in the case V(x) = — aS(x ), the con¬ 
volution is replaced by the sum ( | 20 | . 
even be calculated in O(K) operations. 


Hence, can 


Algorithm 5: calc<i>g 


6. TRIPLE SPLITTING INTEGRATOR 

For our second scheme, we consider the coupled sys¬ 
tem © ,( |18| ) as a whole. From a numerical perspective, 
we have an initial value problem for 


P= 2a/(i?r)E k - 1 k Vk 
k 2 " 

q = 2a/ (Lit) qk 

K— 2 

for k = — k /2 to k /2 — 1 do 

I 7 k(t) = 7fc(0) + r ■ (q k ■ p - Pk ■ q) 


y(t) = 

( ^ | 6 C 2K 

(42) 

7(f) = 

(7_K/ 2 (t) ... 7 «r/ 2 _i(t)) 

(43) 

d(t) = 

(a_jr/ 2 (t) ... aur/ 2 _i(t)) e C K , 

(44) 



(45) 


whose right hand side /(y) can be split into three parts, 
/ (7, a) = Ay + g(a ) + /i(y, a). (46) 


Here, Ay is the first term of the equation of motion 
for a , i.e., 



(47) 


which means that it represents the same action on a as A 
in the nonlinear case above. The function g{a) represents 
the right hand side of the evolution equation for 7 and 
h( 7 , a) is the second term of Eq. (18 1. 


Figure 4: Sketch of the algorithm which for given values 7 ( 0 ), 
<5(0) and a given step size r calculates the solution 7(7") = 
4> Ti 9 ( 7 ( 0 )) to the initial value problem (48 1 . 


6.2. Calculating 


We consider the subproblem 


i = htf(0),a(t)), 

< 3 ( 0 ) = d 0 - 


(50) 


As the right hand side’s Lipschitz constant is small, we 
can apply a standard integration scheme. This yields the 
algorithm outlined in Fig. [5j The appealing fact about 
our triple splitting is that in the case of a contact interac¬ 
tion, the subproblem ( |5C)| ) can be solved exactly in O(k) 
operations as we show now. 

Introducing b € via 

bk = y— (2 7 fe(0) — 1), (51) 
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Algorithm 6: calc<l>/, Algorithm 7: calc h 

Y =calc h( a) c = calc convolution)! 7 , a) 

for k = — k /2 to k /2 — 1 do for k = — K /2 to K /2 — 1 do 


Y(k) = 
a k (t) + r/2Y(k) 


Y =calc_Ji(Y) 
for k = — k /2 to k /2 — 1 do 
a fc (t) = Ofe(t) + r Y (k) 


d = 2 • 7fe — 1 

h(k) = —ic ■ d 


Hence, the solution to the initial value problem (501 is 
given by 


(r) = * T , h (a(0)) = e-^a(0) =: e B a(0). (54) 


We now show that for a given <5(0), c?(t) can be calcu¬ 
lated in G(K) operations. 

For a given n £ N, we have 


Figure 5: The left panel shows the algorithm which for a 
given value <5(nr) and a given time step r calculates a((n + 
1)t) = 4>°™(a(/!T)) with the explicit midpoint rule. The 
right panel shows the algorithm which for a given value a 
calculates h( 7 , Al¬ 


and the K x A-matrix 

b_ K / 2 ■ ■ ■ b_K/ 2 


B = 


we can write 


^ K / 2-1 b K/2-l 


K 




(52) 


(53) 


—i b- K / 2 TC n 1 ... —i b_ K / 2 Tc n 1 

B n = I : '■ : ], (55) 

\—\b K / 2 -iTc n ~ l ... —i6K/ 2 iTc" 1 


K 


with 


: =" iT £ b J 


■? = -TT 


(56) 


Consequently, with Id denoting the A x A identity ma¬ 
trix, 


we have 


exp(—i tB) = Id + ^ — 


/ —i b_ K / 2 Tc n 1 ... —i b_ t 


\-\b K / 2 -\Tc n 1 ... -\b K / 2 -iTc n 1 
( —i6_rr/ 2 r(exp(c) — 1) ... —i6_ K / 2 r(exp(c) — 1) 


= IdH— 

c 


\—i6jf/ 2 _ 1 r(exp(c) - 1) ... -i6*/ 2 _i-r(exp(c) - 1) 


With this, the matrix-vector multiplication in Eq. (541 yields 

( &-^/ 2 (exp(c) - 1) &j{0 ) \ 


exp(— irB)d(0) = <S(0)- 

c 


y&rc/ 2 _i(exp(c) - 1) E / = _4 a j(°)) 


(57) 


(58) 


(59) 


Thus, the solution of the initial value problem (501 can 
efficiently be calculated by the algorithm illustrated in 
Fig. [6] 


Having found efficient algorithms for all three subprob¬ 
lems we have split the system into, we can now recompose 
them. 


6.3. SplitBCS 


As all the three flows § t ,a, and &r,h are at 

least of second order, each symmetric composition of 
them gives rise to a second order integration scheme, see, 




















Algorithm 8: calc<l>^ 

c = ~irJ2Z- 1 i< b * 

K —~2 

s = E E - * «fc(o) 

K — TC 

e = exp(c) — 1 
for k = — K li to k /2 — 1 do 
I Ofc(t) = afc(0) - i r ■ e- s ■ b k /c 


Figure 6: Sketch of the algorithm which for given values 7(0), 
a(0) and a given step size r calculates the solution a(r) = 
47 ,7! (<3(0)) to the initial value problem (501 for the case of a 
contact interaction. 


e.g. ESI Chapter II.5]. We propose the composition 

«J,num = § T AghgA ■- $ T/2 A o <f> T/2 g O <t> T h O <4> T / 2g O 4> t/2jA , 

(60) 


1, have to be calculated. But, if storage is not 
a problem, one only has to calculate the cos and 
sin once at the beginning of the simulation as the 
arguments are the same in each step. This is what 
we did in our implementations. Accordingly, the 
number of operations specified in Tabs. 0 and |TTJ 
refer to this efficient version. 

Putting everything together, we obtain our integrator 

SplitBCS as outlined in Fig. [7] 


Algorithm 9: SplitBCS 


a = 47/ 2 , A (a( 0)) 

for n = 0 to N do 

(7,«) = 37,^9( 7 , a), 
a = $r,A(a). 

a = 4> T / 2 , A (a) 


as this yields the fastest and most accurate scheme among 
the possible combinations as we will see in the next Sec¬ 
tion. If even more accuracy were required, we could use 
a suitable composition of the scheme (601; see [21 125] for 
more information on compositions. 

In the same way as for the algorithm of Section [5j the 
last sub-step of each step can be combined with the first 
sub-step of the following step which reduces the CPU 
effort. Even more computational costs can be saved by 
paying heed to the following points. 


From Eq. (25) it can be deduced that 


d_ 

df 


E *(*)) 

,i=- # / 


= 0 . 


(61) 


Figure 7: Sketch of our algorithm SplitBCS which for given 
initial values 7 ( 0 ), a(0) and a given step size r approximates 
( 7 (Ar),a(Ar)) T = 4>jv-r,/(7(0), a(0)). 


6.4. Number of operations 

In order to compare the efficiency of SplitBCS to the 
one of BCSInt, we count the number of operations re¬ 
quired for the respective sub-algorithms and for SplitBCS 
as a whole, too. The result can be found in Tab. [II] If 
the explicit midpoint rule is replaced by the Cash-Karp 
scheme, the number of operations in the calculation of 
<P T: h increases as for BCSInt. We see that SplitBCS cal- 


Thus, the sum over all 7 k(t), and, as a consequence, 
also the quantities c and e appearing in the calcu¬ 
lation of $ T ,h, cf. Fig. [g] are preserved along evo¬ 
lutions of the equations of motion. Hence, c and e 
only need to be calculated once at the start of the 
simulation when considering a contact interaction. 

• Both and 4 > t ,/i require the computation of the 
convolution, cf. Figs. [4] and [ 6 j However, does 
not modify a which means that the convolution in 
the first call of calc_ft. is the same as the one already 
calculated in calc<l>g. Hence, by suitably combining 
the calculation of 

4* r.ghg ■— 4?-,y 2) g ® ® 

into one algorithm, one can avoid redundancies[27j. 

• The calculation of <1> TjA can be made more efficient 
for both BCSInt and SplitBCS when a fixed step 
size is used. In this case, during each call of <1> r A , 
cos and sin of 2 ( fc2 / i2 — A 4 ) U & = K / 2 ,..., K /2 — 


Algorithm 

^Operations per call 

Calculation of <t> r ,A 

Calc_ h 

calc<£ jj 
calc<fr 9 

Calculation of $> T ,ghg for contact interaction 

SplitBCS 

SplitBCS for contact interaction 

14 • K + 14 

1 x calc_convolution + 7 • K + 12 

2 x calc_ h + 8- K + 9 = 2x calc_convolution + 22 • K + 33 

1 x calc_convolution + 6 • K + 12 

18 • K + 39 

3 x calc_convolution + 34 • K + 53 

32 • K + 53 


Table II: The required number of operations per step as a 


function of dimension of the ODE system (17 1 , (18 1 for the 


sub-algorithms of SplitBCS and for SplitBCS itself. 


culates one convolution more than BCSInt. Thus, BC¬ 
SInt is expected to be faster for general settings with a 
huge number of basis functions. For the important case 
of a contact interaction, however, we are able to calculate 
<!> r.ghg very efficiently. This is why we choose the compo¬ 
sition (601 over other possible sequences of the subflows. 


With this, SplitBCS is even faster than BCSInt in the 
physically important setting. 

Let us now subject the schemes to numerical tests. 
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7. NUMERICAL EXPERIMENTS 

All the numerical experiments presented here were run 
on a Core 2 Duo E6600 machine with 2.4GHz and 4GB 
RAM. In order to have physically realistic data to start 
our experiments with, we chose a system which is slightly 
superconducting. Such a system can be obtained by set¬ 
ting 


tanh ( ^> a -"> 2 +* 3 


2 T 


m _ 1 _ k /L* - H _ 

1 [ 2 2 y/{k 2 /L* - fj,) 2 + h 2 


h 


tanh ^( >!, /y +ha ) 


afe(0) 2 ^(^ - m )2 + h 2 ’ 


(63) 

(64) 


where h. = 0.1 is a small parameter. The critical temper¬ 
ature T of the system depends on the chemical potential 
H and on the interaction potential V. In the simulations 
presented here, we considered a system with a contact 
interaction V[x) = —a8(x). In this case, T can be calcu¬ 
lated from the implicit formula, cf. m, 


2n 

a 


tanh ^ p 2T m 
V 2 ~ P 


dp. 


(65) 


For our simulations, we chose a = fj, = 1 which yields 
T = 0.19. 


As a measure of an integrator’s accuracy, we used the 
discrete energy (27 ) which is conserved along the exact 
solution of the ODE system Thus, the reliabil¬ 

ity of a numerical integration scheme can be checked by 
tracking the relative error A F K , defined by 


A F K {t) 


F K m)M)-F K m),*m 

F K { 7(0), a(0)) 


( 66 ) 


along the numerical evolution. 

We first used this tool to compare SplitBCS to BC- 
Slnt with 4> r j 1 calculated via the fifth order Cash-Karp 
method. For this, we fixed L = 32, K = 256 • L and 
chose a step size r = 0 a /k. We evolved the system until 
t = 0(L ) with both integrators and plotted the rela¬ 
tive error in the energy, A F K , against integration time 
t in the left panel of Fig. [8] We repeated the procedure 
for L = 64 and plotted the result in the right panel of 
Fig-|8l Although the error increases slightly at the end of 
the integration for BCSInt, both schemes seem to be very 
accurate. When comparing BCSInt with < f> r ^ calculated 
via the fifth-order Cash-Karp method to BCSInt where 
& T ,fi w as calculated with the explicit midpoint rule, we 
found no differences in the relative error of the energy. 
Hence, we recommend the use of the latter method as it 
is much faster. 

Other physically relavant constants of motion are the 
eigenvalues A*, of the particle density matrix T. BCSInt 
preserves them by construction. In order to check their 




t 


Figure 8: The relative error A F K of the free energy as a 
function of integration time t for SplitBCS and BCSInt in 
semilogarithmic scale. The left panel shows the result for 
L = 32, the right panel depicts the corresponding result for 
L = 64. 


behavior when using SplitBCS, we also tracked the eigen¬ 
values together with their corresponding relative error, 


A A k(t) 


A k(t) - Afc(O) 
Afc(0) 


(67) 


along the evolution. We found out that, up to very small 
rounding errors, all eigenvalues were preserved for Split¬ 
BCS, too. As an illustration, we plot some eigenvalues 
and the relative error in Ao in Fig. [9] 

With regard to SplitBCS, the question remains as to 
whether we could have done even better by choosing an¬ 
other sequence of the sub-flows than composition (601. 
In order to go into this matter, we also evolved the sys¬ 
tems for L = 32 and L = 64 for various other compo¬ 
sitions of the sub-flows $r,g and <f) T ^h , and again 

plotted A F k as a function of the integration time t. The 
resulting plots are shown in Fig. [TO] We also tested the 
other possible sequences which are not shown in the plots. 
However, we found out that the relative error in the en¬ 
ergy seems only to depend on the spot of 4> T ^ in the 
composition. This means that & T ,AhghA is as accurate 
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t 



t 

Figure 9: The left panel shows some eigenvalues A k of the 
density matrix as a function of integration time t for SplitBCS 
applied to the system with L = 32. The right panel shows the 
corresponding relative error AAo of the density matrix’ first 
eigenvalue in semilogarithmic scale. 




t 

Figure 10: The relative error A F K of the free energy as a 
function of integration time t for SplitBCS and other possible 
compositions in semilogarithmic scale. The left panel shows 
the result for L = 32, the right panel depicts the correspond¬ 
ing result for L = 64. 


as SplitBCS. But we could not find an equally efficient 
implementation for <& Tt hgh as the one for <& T) gh g - This 
is why we strongly recommend the use of the composi¬ 
tion (60 1 , shortly SplitBCS, in simulations of the discrete 
BCS equations with a contact interaction. 

In order to show, as a last point, why standard integra¬ 
tion schemes are of no use for the discrete BCS equations, 
we apply the popular fifth order Cash-Karp scheme of [3] 
to the equations with the same L and the same step size 
as for the splitting methods. When plotting the result- 


as tor the splitting mi 
ing A F k , cf. Fig[lT] 


we observe an exponential growth 
in the error. This is in accordance with theoretical ex¬ 
pectations, see, e.g. m- 

Let us now summarize our results. 


8. CONCLUSION 


In this work, we have presented two fast and accurate 
evolution schemes, BCSInt and SplitBCS, for the coupled 


discrete BCS equations which arise from a Fourier space 
discretization of the BCS equations for superconducting 
materials. BCSInt uses the preservation of the density 
matrix’ eigenvalues to decouple the system and a subse¬ 
quent splitting of the decoupled system into two terms. 
SplitBCS is based on a splitting of the coupled equations 
into three subproblems which for the important case of a 
contact interaction can all be solved exactly by employ¬ 
ing basic operations only. Crucially, the CPU effort for 
these exact solutions grows only linearly in the dimen¬ 
sion of the spatial discretization. Further computational 
costs could be saved by aptly recombining the flows of 
the subproblems. In numerical tests, the schemes have 
been shown to be very accurate. Additionally, they pre¬ 
serve the discrete analog of the physical energy and the 
eigenvalues of the particle density matrix up to very small 
errors. We have, thus, come up with very useful tools for 
simulations in the field of superconductivity. 
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Figure 11: The relative error A F K of the free energy as 
a function of integration time t for the explicit Cash- Karp 
scheme in semilogarithmic scale. The left panel shows the 
result for L = 32, the right panel depicts the corresponding 
result for L = 64. 
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